The standard map

The standard map $f(x,y) = (x+y+a\sin(x),y+a\sin(x))$ is an area-preserving map on the 2-torus $[0,2\pi]^2$ resulting from a symplectic time-discretization of the planar pendulum. For $a = 0.971635$, its phase space shows the characteristic mixture of regular (periodic or quasi-periodic) and chaotic motion. Here, we repeat the experiment in Froyland & Junge (2015) and compute coherent sets:

using Random

a = 0.971635
f(a,x) = (mod(x[1] + x[2] + a*sin(x[1]), 2π),
          mod(x[2] + a*sin(x[1]), 2π))

X = []
for i in 1:50
    Random.seed!(i)
    x = 2π*rand(2)
    for i in 1:500
        x = f(a,x)
        push!(X,x)
    end
end

using Plots
gr(aspect_ratio=1, legend=:none)
scatter([x[1] for x in X], [x[2] for x in X], markersize=1)

Approximating the Dynamical Laplacian by FEM methods is straightforward:

using Arpack, CoherentStructures, Tensors

Df(a,x) = Tensor{2,2}((1.0+a*cos(x[1]), a*cos(x[1]), 1.0, 1.0))

n, ll, ur = 100, [0.0,0.0], [2π,2π]               # grid size, domain corners
ctx, _ = regularTriangularGrid((n,n), ll, ur)
pred(x,y) = peuclidean(x[1], y[1], 2π) < 1e-9 &&
            peuclidean(x[2], y[2], 2π) < 1e-9
bd = boundaryData(ctx, pred)                      # periodic boundary

I = one(Tensor{2,2})                              # identity matrix
Df2(x) = Df(a,f(a,x))⋅Df(a,x)                     # consider 2. iterate
cg(x) = 0.5*(I + dott(inv(Df2(x))))               # avg. inv. Cauchy-Green tensor

K = assembleStiffnessMatrix(ctx, cg, bdata=bd)
M = assembleMassMatrix(ctx, bdata=bd)
λ, v = eigs(K, M, which=:SM)

using Printf
title = [ @sprintf("\\lambda = %.3f",λ[i]) for i = 1:4 ]
p = [ plot_u(ctx, v[:,i], bdata=bd, title=title[i],
             clim=(-0.25,0.25), cb=false) for i in 1:4 ]
plot(p...)